function matrix_A = assemble_matrix_A_2d(node, elem, function_cor, basis_type, trial_type, test_type)

gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
gauss_weights_rec = [0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.79012345679012345679012345679012, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864];

matrix_A = sparse(size(node,1), size(node,1));

phi_x = reference_basis_points_2d(basis_type, "dx");
phi_y = reference_basis_points_2d(basis_type, "dy");
psi_x = reference_basis_points_2d(basis_type, "dx");
psi_y = reference_basis_points_2d(basis_type, "dy");
ref_trans_tri = reference_basis_points_2d("P1");
ref_trans_rec = reference_basis_points_2d("Q1");

switch basis_type
    case {"P1", "P1b", "P2"}
        for n = 1:size(elem,1)
            E = node(elem(n,:),:);
            V = E(1:3,1:2);
            J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
            switch trial_type
                case "dx"
                    phi = cellfun(@(x,y) (V(3,2)-V(1,2))/J_det .* x + (V(1,2)-V(2,2))/J_det .* y, phi_x, phi_y, "UniformOutput", false);
                case "dy"
                    phi = cellfun(@(x,y) (V(1,1)-V(3,1))/J_det .* x + (V(2,1)-V(1,1))/J_det .* y, phi_x, phi_y, "UniformOutput", false);
            end
            switch test_type
                case "dx"
                    psi = cellfun(@(x,y) (V(3,2)-V(1,2))/J_det .* x + (V(1,2)-V(2,2))/J_det .* y, psi_x, psi_y, "UniformOutput", false);
                case "dy"
                    psi = cellfun(@(x,y) (V(1,1)-V(3,1))/J_det .* x + (V(2,1)-V(1,1))/J_det .* y, psi_x, psi_y, "UniformOutput", false);
            end
            ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
            for alpha = 1:size(elem,2)
                for beta = 1:size(elem,2)
                    f = (function_cor(ref_points(:,1), ref_points(:,1))' .* phi{alpha} .* psi{beta});
                    r = (gauss_weights_tri * f') .* abs(J_det);
                    matrix_A(elem(n,alpha), elem(n,beta)) = matrix_A(elem(n,alpha), elem(n,beta)) + r;
                end
            end
        end
    case {"Q1", "Q1b", "Q2"}
        for n = 1:size(elem,1)
            E = node(elem(n,:),:);
            V = E(1:4,1:2);
        end
    otherwise
        error("Invalid basis type");
end

end

